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ABSTRACT 

We present dual-wavelength observations and modeling of the nearly edge-on 
Class young stellar object LI 157- mm. Using the Combined Array for Research 
in Millimeter-wave Astronomy, a nearly spherical structure is seen from the cir- 
cumstellar envelope at the size scale of 10^ to 10^ AU in both 1 mm and 3 
mm dust emission. Radiative transfer modeling is performed to compare data 
with theoretical envelope models, including a power-law envelope model and the 
Terebey-Shu-Cassen model. Bayesian inference is applied for parameter estima- 
tion and information criteria is used for model selection. The results prefer the 
power-law envelope model against the Terebey-Shu-Cassen model. In particular, 
for the power-law envelope model, a steep density profile with an index of ~2 is 
inferred. Moreover, the dust opacity spectral index /3 is estimated to be ~0.9, 
implying that grain growth has started at LI 157- mm. Also, the unresolved disk 
component is constrained to be < 40 AU in radius and < 4-25 Mj^p in mass. 
However, the estimate of the embedded disk component relies on the assumed 
envelope model. 
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Introduction 



Protostars are surrounded by their natal envelopes in the earliest stage of evolution. 
These envelopes supply the material that is actively infalling onto the embedded star-disk 
system, and their properties can affect the subsequent evolution. Much theoretical work 
has been done to address the collapse process and the physical properties of the envelope, 



netic fields (e.g., 


Larson 


1969; 


Penston 


1969; 


luntei 


1977; 


Shu 


1977 


Terebev et al. 


1984 




Whit 


worth & Summers! 


1985; 


Galli & Shi. 


1993 


Masunaea & Inutsuka 


2000; 


Tassis & Mouschovias 


2005; 


Hennebelle & Fromane 


2008 


). While various theoretical models of the collapsing en- 



velope have been suggested, distinguishing between the suggeste d theoretical models has 
alway s been an observational challenge. For example, the model of iTerebey. Shu, fc Cassen 
( 1 1984 hereafter the TSC model) has been widely used and consistent results have been 
obtained, esp ecially in fitting the s pectral energy distributions of unresolved young stellar 



objects (e.g., iRobitaille et al.l 120071 ). but it is based on the Shu (1977) model, which could 



not fit a sample of Class protostars with reasonable ages. (Looney et al. 2003). 

To better constrain the envelope structure, we carry out complete modeling with dual- 
wavelength millimeter data. At such wavelengths, the continuum is dominated by dust 
emission from the envelope and the embedded disk. Interferometry is a useful tool to probe 
the structure of protostellar envelopes as it measures emission at various spatial scales, 
leading to a more complete analysis for the envelope. By comparing the predicted envelope 

"is, theoretical collapse models can be tested (e.g.. 



Chiane et al. 


2008; 


Maurv et al. 


2010) 



us to constrain the physical properties of the embedded disk component. 

In this paper, we focus on the edge-on Class protostar L1157-mm (also k nown as LI 157- 



1992; 


Kun 


1998; 


Kun et al. 


2008) 



A chemically active outflow driven by L1157-mm has be en detected in multiple species 



( iGueth et al.l Il996l ; iBachiller et al.l l2001t iNisini et al.l 120101 ). Perpendicular to the outflow 



orientation, a flattened envelope with a linear size of ~20,000 AU is seen in 8 nm absorption, 
N2H+ emission, and NH q emission, showing complex kinematics from rotation, in fall, and 



outflow in the envelope (iLooney et al.l 120071 : IChiang et al.ll2010l ; iTobin et al.ll201l[ ). On the 



other hand, dust contin uum traces the envelope structure as well as reveals a compact core 
(e.g.. lGueth et al.ll2003l ). Additionally, the presence of a circumstellar disk embedded inside 



the e nvelope is suggested by methanol observations (iGoldsmith et al.lll999l ; IVelusamy et al. 
2002h . 



We have collected 1 mm and 3 mm interferometric data at multiple array configu- 
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rations using the C ombined Array for Research in Milhmeter-wave Astronomy (CARMA; 
Woody et all booi l @. g2] gives an overview of the observations and the data reduction. 



Details of the modehng are addressed in §3] with a further supplement on our statistical 
approach in the Appendix. The results are presented in §11 their implications are discussed 
in ^ and a summary is given in ^ 



2. Observations and Data Reduction 

L1157-mm was observed by the 15-element CARMA between Oct 2007 and Jan 2010. 
At that time, the science array of CARMA consisted of six 10.4-meter antennas and nine 
6.1-meter antennas. Dust continuum at both 1 mm and 3 mm bands was observed using 
multiple array configurations, as summarized in Table [H 

The phase center of the observations before September 2008 was a = 20'^39™06';20, 5 = 
68°02'15'.'9, and shifted to a = 20''39™06!26, 6 = 68°02'15'.'8 afterwards as more precise coor- 
dinates were determined by high resolution observations. However, all data presented here 
have been corrected to have the common phase center at a = 20^39™06!26, 6 = 68°02'15'.'8 
(J2000). 

The main phase calibrator for all tracks was 1927+739 (with the exception of one A- 
array track) and was observed with a phase calibrator-source cycle of 10-15 minutes. For 
all A- and B-array observations, a weaker quasar, 2009+724, was observed as the secondary 
phase calibrator. The secondary phase calibrator was not used in the calibration process; 
instead, it is used to check the point source response. For 3 mm A- and B-array tracks 
observed in winte r 2009-2010, the CARMA Paired Antenna Calibration System (C-PACS; 



Perez et al.ll2010l ) was employed to calibrate the atmospheric phase variation on short time- 
scale. With C-PACS, a reference array continuously monitors a nearby quasar, called the 
atmospheric calibrator, for atmospheric delay, while the science array observes the science 
target. The eight 3.5-meter antennas, from the previous Sunyaev-Zel'dovich Array (SZA), 
were used as the reference array. The C-PACS correction is effective for data with long 
baseline. 

The data reduction, calibration, and imaging were done usin g the Multichanne l Image 
Reconstruction, Image Analysis and Display package (MIRIAD; Sault et al. 1995 )E The 



bandpass and flux calibrators for each track are listed in Table [H After the data are reduced. 



^http://www.mniarray.org/ 
■^http://carnia. astro. umd.edu/miriad/ 
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Table 1. Summary of Observations 



Frequency 


Array 


Date 


Observing 


Bandpass 


Flux 


Beam Size ^ 


Beam P.A. 


(GHz) 


Config. 




Time (hr) 


Calibrator 


Calibrator 


(") 


(degree) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


229 


13 


2007-12-17 


1.5 


/" ^ ATA 

3C454.3 


TV 4"\ li T/^ A f\ 

MWC349 


0.4x0.3 


-75 




c 


2008-04-13 


3.1 


3C454.3 


MWC349 


1.0x0.8 


-66 




D 


2008-03-07 


3.1 


3C454.3 


TV /TTir/^O Ad 

MWC349 


2.3x2.0 


-29 


91 


A 

A 


2009-01-27 


3.3 


1642-1-689 


IV /fViV r/^o Art 

MWC349 


0.4x0.3 


-65 






2010-01-26 


2.2 


3C273 


MWC349 










2010-02-01 ^'^ 


2.9 


3C454.3 


Neptune 








B 


2007-11-17^ 


4.5 


175H-096 


MWC349 


0.9x0.8 


-83 






2007-11-19 ^ 


2.4 


3C454.3 


Neptune 










2007-11-20 ^ 


1.5 


3C273 


3C273 










2009-12-14 b 


0.6 


3C454.3 


Uranus 










2009-12-15 '''^ 


5.4 


3C345 


MWC349 








C 


2007-10-03 


2.5 


1751-h096 


MWC349 


2.3x2.0 


-87 






2007-10-05 


2.8 


1751-1-096 


MWC349 










2008-04-05 


4.8 


3C273 


MWC349 








D 


2008-02-29 


3.1 


3C454.3 


Uranus 


6.0x5.2 


84 






2009-03-19 


6.0 


3C345 


MWC349 










2009-03-20 


6.2 


3C345 


MWC349 










2009-03-27 


0.8 


1642+689 


MWC349 










2009-03-29 


3.0 


1642+689 


MWC349 








E 


2008-10-02 


4.2 


3C454.3 


Uranus 


11.5x10.2 


78 






2008-10-05 


3.5 


3C84 


MWC349 







"^The synthesized beam of the combined data with natural weighting at each array configuration 

'"Track examined with the secondary phase calibrator 2009+724 

'^Track using the primary phase calibrator 1849+670 instead of 1927+739. 

•^Track calibrated with C-PACS using the atmospheric calibrator 2022+616 
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the flux density of botli primary and secondary pliase calibrators is plotted as a function of 
u-v distance in order to verify a flat trend, implying that decorrelation is not significant at 
long baselines. 

The largest uncertainty of interferometric data comes from flux or absolute amplitude 
calibration. Although independent of relative brightness and image quality for one track of 
data, it affects the analysis through the differences between tracks. The uncertainty can be 
larger than 10% mostl y due to the planetary mod el of the flux calibrator used in the data 



reduction process (e.g.. lMoreno fc Guilloteaull2002l ). The large uncertainty cannot be avoided 



unless the planet modeling is improved. The flux uncertainty can affect the analysis through 
(1) the uncertainty between tracks at the same wavelength, and (2) the uncertainty between 
tracks at different wavelengths. The uncertainty of the first kind can affect the deduced 
envelope and disk structure in the modeling. To ensure its impact is minimized, we compare 
the flux of the common phase calibrator 1927+739 among tracks. At each wavelength, we 
verify that the flux value varies smoothly in time and is consistent with the flux reported 
in the standard CARMA/MIRIAD catalog. Also, data and model are compared at each 
visibility point, so the uncertainties are better preserved compared to modeling using binned 
visibilities (see §3.5p . For the rest of the paper we consider the absolute flux uncertainty of 
the second kind for dual-wavelength analysis. A 10% uncertainty for the absolute flux at 
each wavelength is adopted. It dominates errors in estimating spectral parameter, such as 
the dust opacity spectal index, as will be seen in §11 Other model parameters can be affected 
directly or indirectly through the uncertainty of the spectral parameter. Other systematic 
uncertainty from instruments and calibrations may exist and propagate in the analysis as 
well, but are presumably less than 10%. 

The reduced data of the science target L1157-mm are shown in Figured] by the annuli- 
averaged flux density with respect to u-v distance. Continuum data from all spectral windows 
are combined. Figure [2] presents the continuum maps of L1157-mm. Super-uniform weight- 
ings with different robustness parameters are used to obtain different synthesized beamsizes 
in order to emphasize envelope structures at different size scales. The continuum of LI 157- 
mm shows spherical structures from 2000 AU scale (8") down to 100 AU scale (0.4"). In panel 
(f), the envelope structure is slightly elongated perpendicular to the ou tflow direction, but 



the larger-scale extended envelope, detected by IRAM 30-m telescope (iGueth et al.ll2003l ) 
and SMA (with lower resolution than our CARMA observations; Tobin et al. in preparation), 
is not seen in our CARMA dust continuum data. No apparent disk or flattened structure is 
seen at small scale either. 
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Fig. 1. — Flux density of L1157-mm dust continuum at 1 mm {filled symbols) and 3 mm 
{open symbols). The visibilities are vector averaged around the source center and binned into 
u-v annuli. Data collected in different array configurations are plotted separately. The error 
bars show only the statistical errors within each annuli-bin, while the typical uncertainty 
carried by single data visibility is around 0.2-2.5 Jy. 
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Right Ascension offset (arcsec) 

Fig. 2. — CARMA 3 mm (upper) and 1 mm (lower) dust continuum images of LI 157- 
mm. The same multi-configuration data with different u-v imaging weightings are shown to 
emphasize structures on different size scales. The contour levels, noise rms (a), and beams 
are as follows: (a) [3,4,5,7,10,14,18,22] xcr, a = 0.9 mJy beam-^ 2.40"x2.03" at a position 
angle of 90°; (b) [3,4,5,7,10,14,18,22,26] xcr, a = 0.6 mJy beam-\ 1.34"xl.l0" at a position 
angle of -88°; (c) [3,4,5,7,10,14,18] xo", a = 0.6 mJy beam-^ 0.65"x0.54" at a position angle 
of -82°; (d) [3,4,5,6,7] xcT, a = 0.9 mJy beam-\ 0.32"x0.28" at a position angle of -73°. 

(e) [3,4,5,7,10,14,20,30,42] XCT, a =4.0 mJy beam^^ 1.77"xl.61" at a position angle of -38°; 

(f) [3,4,5,7,10,14,18,22,26] xa, a =5.5 mJy heam-\ 1.14"x0.98" at a position angle of -48°; 

(g) [3,4,5,7,10,13,16] XCT, a = 7.0 mJy beam"!, 0.71"x0.62" at a position angle of -51°; (h) 
[3,4,5,6,7,8] X(T, a = 12.0 mJy beam-^ 0.37"x0.30" at a position angle of -83°. 
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3. Aspects of Modeling a Class YSO 

To compare a YSO model with observations, we consider the physical conditions of the 
system, including the density and temperature structures ( §3.11 and §3.2p and dust grain 
properties ( §3.3p . The radiative transfer tool RADMC-3D, developed by C. P. DuUemond 
and co-authors ( DuUemond &: Dominik 20041 ) @, is used. Observational effects from the 



interferometer are taken into account and a Bayesian approach is taken for model fitting 
( §3.5p . In the following sections, we discuss the details of each facet in the modeling. 



3.1. Envelope Structure 



A simple Class YSO model consisting of a spherical dusty envelope, a bipolar outflow, 
and possibly a circumstellar disk is considered. For the envelope structures, we examine a 
simple power-law density profile, representing self-similar collapse solutions, and a collapse 
with rotation (Terebey et al. 1984). An unresolved component is included to represent a 
compact disk structure. 

To include a simple bipolar outflow cavity in the model, we remove material orientated 
with the observed envelope geometry and outflow properties: a position angle of 152° for the 
outflow-axis c ut, an inclination angle of 80°, and an ope ning angle 30° for the outflow cavity 
are assumed (jChoi et al.lll999l : iGueth et al.lll996l . 119971 ). For simplicity, the inner and outer 
radii of the envelope are fixed to be 12 AU and 10,000 AU, respectively. The inner envelope 
cavity is smaller than the highest observational resolution, and always within the central cell 
in the model images. A large outer radius is adopted so there is no ringing effect due to 
interferometric response on a sharp cutoff in the envelope. Additionally, as the density and 
temperature are much lower in the outer envelope, precise choice of the outer radius does 
not play an important role at these wavelengths. 



3.2. Temperature Structure 



While many theoretical models ignore the internal heating from the newborn protostar, 
it is critical to take into accoun t the protostellar contribution to agree with observational 
luminosity (lAdams fc Shulll985l ). The heating and cooling of dust grains, dominated by 
the central illumination and dust grain properties, should be balanced to obtain an equi- 



"'http: / / www.ita.uni-heidelberg.de/~dulleniond/software/radnic-3d/ 
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librium temperature. To simulate millimeter-wave observations of protostars surrounded by 
dusty environments, such a realistic temperature distribution needs to be either assumed or 
calculated. 

The temperature structure can be approximated assuming simple conditions of the dusty 
envelope. Assuming a centrally illuminated spherical envelope in which the density has a 
power-law dependence on radius, 



where po is the density at an arbitrary radius tq and p is the density power-law index, and 
assuming a pure power-law dust opacity with a spectral index /3, the temperature structure 
in the optically thin outer envelope can be approximated by 



T(r) 



2 



(2) 



(jWolfire fc Cassinellilll986l : lAdamslll99ll ). 



However, if the assumptions of power-law dust opacity and spherical power-law density 
do not hold, the approximation in Eq. ([2]) can be inadequate even in the optically thin reion. 
For example, the dust opacity is not a pure power-law at short wavelengths, and the density 
structure can also be more complicated than the power-law profile. Furthermore, Eq. ([2]) is 
only valid in the optically thin region and relies on Tq at tq. With a fixed central heating 
source, Tq at tq is characterized by the optically thick-thin transition zone, and is difficult 
to estimate without a good understanding of the optically thick inner region. 

Given the difficulty to approximate the temperature structure with a variety of envelope 
models and ranges of model parameters, we calculate a self-consistent temperature distribu- 
tion for each set of paramete rs using the Monte Carlo radiative t ransfer code RA DMC-3D 
( iDuUemond fc Dominikll2004j ). A luminosity of 8.4 Lq is adopted (jFroebriclj|2005l ) as a fixed 
input in the radiative transfer calculation. This bolometric luminosity can be underestimated 
due to insufficient sampling of the spectral energy distribution, but can be overestimated 
due to a larger assumed distance. Furthermore, the intrinsic luminosity can be larger than 
the measured bolometric lurn i nosity due to the sou rce's edge-on orientation (e.g., see more 
discussions in lFroebrichll2005l : IWhitney et al.l 120031 ). Despite the uncertainty in luminosity, 
a self-consistent temperature structure is the best compromise for now. 
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3.3. Dust Grain Properties 



Dust grain properties such as chemical composition, geometry, alignment, degree of ion- 
ization, and size distribution play important roles in star-forming processes from thermody- 
namics and grain surface chemistry to timescales of magnetic fi eld effects. For dust grains in 
the di ffuse interstellar medium, the classic model constructed bvlMathis. Rumpl. fc Nordsieck 
( 1977 . hereafter MRN) with optical constants calculated by Draine fc Leel ( 1984 ) can repro- 
duce the interstellar extinction and polarization observations from infrared to ultraviolet 
wavelengths. However, for dust grains in dense cores and star forming regions, collisions and 
interactions between grain particles become more important. lOssenkopf fc HenningI (119941 ) 
has considered the dust coagulation process in dense protostellar cores and found that the 
opacity can be enhanced by a factor of a few as grains aggregate. The authors started with 
the MRN grains covered with different amounts of ice mantles, and investigated the optical 
constants after 10^ yrs of coagulation in gas densities ranging from 10^ to 10^ cm~^. 

In our modeling, we adopt the dust opacity, or the mass a bsorption coefficient k define d 
as the cross section per unit mass, from column 5 of Table 1 in lOssenkopf fc Henning] ( 1l994j ). 
the so-called 0H5 grain which is covered by a thin layer of ice mantle and coagulated at 10^ 
cm~^. Besides being widely used, the 0H5 model shows agreements with, and in some cases 



1999; 


Evans et al. 


2001; 


Shirlev et al. 


2005, 


2011a) 



13 



(3) 



At far-infrared and millimeter wavelengths, k can be approximated as a power law with 
respect to frequency 

K= [ — 

This sub-millimeter dust opacity spectral index /3, which can only be studied with multi- 
wavelength observations, varies with environment and is related to grain properties men- 
tioned previously. (3 is > 1.7 in the diffuse interstellar medium and starless cores (e.g.. 



Beckwith fc Sargent 



1991 



Draine fc Leelll984j;ISchnee et al.ll2010[). but signifi c antly lower in proto planetary disks (/? < 1 



Natta et al.l 120071 : iRicci et al.l l2010bl ) . One explanation for 



e.g., 

lower /3 is a larger grain size and more discussions will be in §5.4. In order to better un- 
derstand the dust property of L1157-mm, we include /3 as a model parameter. Based on 
the 0H5 model, we modify the opacity curve with a power-law of index f3 as in Eq. ([3]) at 
wavelengths longer than an arbitrary choice of 700 fim. 

The dust opacity used in the analysis can substantially affect the deduced spectral 
energy distribution as well as temperature structure of YSOs. With the radiative transfer 
tool, we find that the temperature structure is mostly determined by the dust opacity at 
short wavelengths. /3, which characterizes the dust property at long wavelengths, plays a less 
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important role for the temperature structu re; instead, /3 affects the observed flux directly 
through dust opacity (IChandler et al.lll998[ ). 



With the optically thin assumption, /3 can be estimated using flux ratio between two 
wavelengths in the Rayleigh- Jeans regime; here we call the approximate dust opacity spectral 
index /3thin- In this limit, the flux density F^, oc KuBi, oc z/^*'^'"+^; therefore. 



thin 



InFi -InFa 
In Ui — In z/2 



(4) 



[e.g., 



Kwon et al.ll2009[ l. Figure [3] shows Pthin of LI 157- mm using the annuli-averaged vis- 
ibility at each u-v distance bin of our dual-wavelength data. The uncertainty of the (3thin 
estimation is discussed in Appendix A. (3thin is a good approximation in t he optically thin 



region, but a correction te rm is needed in the optically thick region (e.g., iRodmann et al. 
20061 ; iLommen et al.l 120071 ). To avoid the need of the correction term, full optical depth ef- 
fect is considered in our radiative transfer modeling. Nonetheless, Pthm provides a quick and 
rough estimate for the (3 value across the envelope and reveals possible radial dependence. 
As seen in Figure [3l no strong radial dependence is suggested for /3 at LI 157- mm. Therefore, 
we assume uniform grain properties in the envelope model for simplicity. In other words, k 
is only a function of frequency and independent of radius in our envelope model. 



3.4. Free-free Contamination 



We ignore the contribution of free-free emission in this study. Free-free emission from 
ionized winds or jets can contribute partial flux at millimeter wavelengths (~20% at 7 mm, 
Rodmann et al.ll2006l ) and affect model parameter estimates especially for (3 and disk com- 
ponent. However, it plays a minimal role for our data of L1157-mm at 1 mm and 3 mm. 
By extrapolating fluxes at 8.5 GHz and 4.86 GHz (iMeehan et al.lll998l ) to our observed fre- 
quency, we estimated the free- free emission to be around 0.53 mJy at 3 mm and 0.39 mJy 
at 1 mm for L1157-mm. The free-free correction is negligible in the analysis. 



3.5. Model Fitting 



Given a set of model parameters, we estimate the sky brightness distribution for the 
dust continuum with radiative transfer calculations. The density and self-consistent temper- 
ature distributions in three dimensions are considered along with the model dust property. 
Essentially, for each pixel on the plan e of sky, the fl ux is calculated by integrating the dust 



emission along the line of sight (e.g., lAdamsl 1 19911 ). Assuming no background brightness. 



2500 AU 
10" 




1.2 1.4 1.6 1.8 2 

Iog10(uv distance (k^)) 



2.2 2.4 2.6 2.8 



Fig. 3. — Approximate dust opacity spectral index l^thin of L1157-mm as a function of 
u-v distance assuming the optically thin condition. The error bars show only the statistical 
errors without the absolute flux uncertainty, while the shade shows the errors including the 
absolute flux uncertainty. 
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the specific intensity can be expressed as 

I,= [ B,{T) e-^^ dr, = f B,{T{r)) e'^"^'^ p{r) K,dr, (5) 

J los J los 

where By{T) is the Planck function at dust temperature T, p is the envelope density, r 
denotes the position, and Ty is the optical depth from the position r along the line of sight 
{los) to the observer 

Ti,{r) = Ki, I p[f)df= / p{r)dl'. (6) 

Jlos Ji 

T, p, and are all dependent of r. 

With the model sky brightness, we simulate interferometric observations and generate 
model visibilities. The sky image is convolved with the primary beam patterns of the anten- 
nas and then Fourier transformed into visibilities with the observational u-v sampling. In 
the case of the 15-element CARMA, the 6.1-meter dishes and 10. 4- meter dishes give 3 types 
of baselines. Therefore we construct separate primary-beam-corrected images for each kind 
of baseline, and sample the images with corresponding u-v spacing for each data visibility 
from real observations. In addition, images at two wavelengths are constructed individually 
based on the same model. 

Model visibilities are compared with observational data at each u-v sample and wave- 
length. The analysis is done in the visibility domain so as to avoid the complexity brought 
by the CLEAN algorithm, u-v sampling, and imaging process. Some information is lost in 
the image domain through the imaging process, since structures in images can be sensitive 
to beamsize or weighting. In other words, emission at different size scales can either be em- 
phasized or suppressed, causing biases in the model-data comparison; therefore, we perform 
the analysis in the visibility domain. Furthermore, visibilities are com pared data point by 



data point; no binning nor averaging are done (e.g., Ilsella et al.ll2009[ ). 



Assuming the noise from observations is normally distributed or Gaussian noise, the 
goodness of a model-fit can be characterized by the standard chi-square statistics. Real and 
imaginary parts of each visibility point are considered individually, as in 

~{Re{Vmodel,i) - Re{Vdata,i)y ^ {Ifn{Vmodel,i) " Im{Vdata,i)y 

a, 



2 l-n-ei^ Vmodei,i; — -n-ei^ Vdata.ij; -T yil I lyv model,i — -L' ' l-y^ data,i 

^ =2^ ^ — -2 (7) 



where i stands for each visibility point at its unique u-v. The noise of each visibility a is 
the square root of data variance (outputted by MIRIAD task uvinfo) multiplied by a scaling 
factor to account for imperfect weather conditions. The noise level before the scaling follows 



On 



int 
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where kj, is the Boltzmann constant, Tsys is the system temperature, rja is the aperture 
efficiency, rj^ is the correlator efficiency, A is the antenna collecting area, is the number of 
antennas, Au is the bandwidth, and tint is the on-source integration time. 

The scaling factor is used to correct a for the phase decorrelation and is 1 if no scaling is 
done. There are multiple ways to scale the noise, and the scaling factor should be somewhat 
dependent on the baseline length. In this work we determine the scaling factor by the phase 
scatter in each array configuration at each wavelength. Nonetheless, the factor is always 
larger than 1; in other words, we only adjust cr to make data less constraining. 

W e take the Bayesian approach to compare data and model fe.g.. lFordl[2005[ ISpergel et al 



20071 ). Given a specific model, a global minimum of is searched and verified to be a good 
fit with a chi-square hypothesis test. Then, the Markov chain Monte Carlo (MCMC) method 
is utilized to calculate the posterior probability distributions; in particular, the posterior- 
weighted value and the uncertainty are estimated for each model parameter. Details of our 
ffiting technique are discussed in Appendix B. Note that the deduced parameters are valid 
only within the framework of model assumptions. Evaluating the goodness of a model and 
comparisons between models will be discussed in §5.11 



4. Results 



In this section, the modeling results based on the details described in ^ are presented. 
Three models are considered and shown individually. 



4.1. Spherical Power-law Envelope Model 

We first consider a spherical envelope with a power-law density profile and self-consistent 
temperature structure. In this simplest model, three model parameters are included: (a) the 
dust opacity spectral index /3 as in Eq. ([3]), (b) the dust density po at 100 AU, which scales 
with the total envelope mass, and (c) the density power-law index p as in Eq. ([1]). All other 
model properties are fixed as described in ^ 

We begin with only considering the statistical noise of data visibility and ignoring the 
uncertainty of absolute flux. This is to characterize the case without absolute flux uncer- 
tainty, as well as demonstrate the effect of absolute flux uncertainty. Using the MCMC 
results, the expectation values and uncertainties of all parameters are estimated (Table 
and the marginalized posterior probability distributions in 1-D and 2-D parameter space are 
shown in Figure |H We choose to list the radius of the 68% confidence interval in Table [2] as 
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it represents 1 a. The 68% and 95% confidence regions are shown by 2-D contours to reveal 
any correlation between parameters. For example, as shown by the marginalized contours in 
the 13- Pq plane in Figure S] (middle row, left column), a correlation between the dust opacity 
spectral index /3 and the density scaling po is implied. This is because a larger /3 means a 
smaller dust opacity k in the model, resulting in larger deduced mass and thus larger density 
scaling. 

The parameters are determined with a high precision within the framework of model 
assumptions. The narrow uncertainties can be understood since approximately five million 
independent data visibilities are used to fit only three model parameters. Strong assumptions 
are impo sed in the model. Similar results have been obtained in other studies as well. For 



example, iKwon et al.l (120111 ) have taken the Bayesian approach to estimate model parameters 
and their errors in the applications of T-Tauri disks, and small uncertainties are obtained 
when only the statistical errors in the data are considered. 

However, the absolute amplitude uncertainty, originated by the absolute flux calibration 
in the data reduction process, brings more uncertainties to the model parameter estimation. 
As we have verified that the calibrator flux is consistent among all observational tracks at 
multiple array configurations O, the absolute flux errors effectively cause an uncertain 
scaling to the amplitude of all data at each wavelength. The absolute flux uncertainty can 
play a dominating role in estimating frequency-dependent parameters, but cause minimal 
effects at the relative spatial structures probed at one single wavelength. 

Marginalization in the framework of Bayesian statistics allows us to quantitatively take 
the absolute flux uncertainty into consideration. We introduce two additional nuisance pa- 
rameters, Simm and Szmm, to scale the absolute amplitude of all data at 1 mm and 3 mm, 
respectively. Simm = 1 and S^rnm = 1 means no scaling is done, as in the presented dataset. 



Table 2. Power-law Envelope Model 



Parameter 


Mean 


Radius of the 68% 


confidence interval 






(statistical noise only) 


(with flux uncertainty) 


(1) 


(2) 


(3) 


(4) 


dust opacity spectral index /3 


0.84037 


0.00014 


0.11 


density at 100 AU (in 10^^*^ g cm~^) 


8.588P 


0.0013 


1.2 


density power-law index p 


2.00939 


0.00020 


0.03 



corresponding to a total envelope mass of 1.8120 
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Fig. 4. — Marginalized posterior probability distributions for the parameters of the power- 
law envelope model (dust opacity spectral index (3, density po at 100 AU in the unit of 
10~^^ g cm~^, and envelope density power-law index p). Only the statistical errors, but not 
the absolute flux uncertainty, are considered for the data. In the histograms, the dashed 
vertical lines enclose 68% or 1 a confidence interval, with the expectation values and a listed 
in Tabled The dark and light areas in the 2-D contour plots are the 68% and 95% confidence 
regions. 
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Since the main uncertainty of flux calibration results from the choice of planetary models 
and no model is preferred, a flat probability distribution for both Simm and S^rnm, ranging 
from 0.9 to 1.1, is assumed. The range of the scaling factors is chosen to be consistent with 
the commonly q uoted 10% e rrors for the absolute flux calibration. Our approach is similar 



to the method of iLay et al.l (119951 ). 



Figure [5] shows the marginalized posterior probability distributions of model parameters 
with consideration of the absolute flux uncertainty. The uncertainties of parameters are listed 
in Table [2] column 4. Inclusion of absolute flux uncertainty increases the parameter errors 
by a factor of 2-3 orders of magnitude, and it is critical for parameter estimation as it makes 
data much less constraining. In particular, the dust spectral index /3 is mostly determined 
by the flux ratio between 1 mm and 3 mm, hence it becomes much more uncertain due to 
the uncertainty of absolute amplitude. 

For additional visualization. Figure [6] shows the observational data visibility of L1157- 
mm and the model calculated with the marginalized parameters as an a posteriori com- 
parison. Because there are about five million data visibilities and each visibility contains 
low signal-to-noise, plotting them all does not show information. Therefore, visibilities are 
averaged vectorially and binned in u-v annuli around the source center using the MIRIAD 
task uvamp. The error bars in Figure [6] are statistical errors in the bins. Note that visibility 
data are not averaged in the modeling process, and the binned visibility is just one data 
representation. The same dataset can have multiple representations depending on how they 
are binned, and the statistical errors in the bins may not reflect the whole uncertainty. 

Figure [7] and Figure [8] continue the a posteriori check and compare the model with 
the data in the image domain for the 3 mm and 1 mm dust continuum, respectively. We 
image the model visibilities in the same way the data visibilities are imaged as shown in 
Figure [21 that is, the same sets of u-v imaging weightings are used for showing structures 
at four size scales at each wavelength. Residuals in the visibility domain are also imaged 
and shown in Figure [7] and Figure [8] to demonstrate the fitting error in the image space. 
The subtraction of model from data leaves no residuals greater than 3 a level at the 3 mm 
images, confirming that a good fit is obtained. In the large-scale image of 1 mm continuum, 
the residuals extend towards the north-west of the protostar, which aligns with the outflow 
direction and is likely due to the asymmetric structure from the outflow. In the small-scale 
image of 1 mm continuum, a 3 o" peak is seen to the north-east of the protostar, which is 
likely caused by the differences of the emission peak position measured using 1 mm and 3 
mm data. We estimate the protostar position by fitting a Gaussian to the highest resolution 
observations at 3 mm, and there is a slight offset relative to the protostar position measured 
using 1 mm data. If we shift the model with this offset at 1 mm, no residuals higher than 3 
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Fig. 5. — Same as Figure H] but with the consideration of the absolute flux uncertainty, a 
for the model parameters are listed in Table [2] column 4. 



- 19 - 



o 
o 
o 



>^ 

E o 



(0 

c 

Q 
X 



100 
— P"" 



(arcsec) 
10 

I 



1mm 



® 







3mm 




10 100 
uv distance (kA) 

Fig. 6. — Flux density of the observational data (circles for 1 mm data and asterisks for 3 
mm data) and the model fit (solid lines) with the marginalized parameters for the power- 
law envelope model. While modeling is done with non-averaged visibilities, annuli-averaged 
visibilities are shown as a function of u-v distance. Error bars are statistical errors in the 
annuli-bins only and different from the visibility uncertainty used in the modeling. 
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Fig. 7. — Comparison between 3 mm dust continuum data {upper row, as shown in Figure 
|2]), model {middle row), and residuals {lower row) of L1157-mm in the image space. The 
power-law envelope model is used. Images in each column share the same u-v imaging 
weighting and contour levels. The contour levels, noise rms (a), and beams are: column 
1: [-3,3,4,5,7,10,14,18,22] X(T, a = 0.9 mJy beam-\ 2.40"x2.03" at a position angle of 90°; 
column 2: [-3,3,4,5,7,10,14,18,22,26] xa, a = 0.6 mJy beam-\ 1.34"xl.l0" at a position 
angle of -88°; column 3: [-3,3,4,5,7,10,14,18] xa, a = 0.6 mJy beam-\ 0.65"x0.54" at a 
position angle of -82°; column 4: [-3,3,4,5,6,7] x a, a = 0.9 mJy beam"\ 0.32"x0.28" at a 
position angle of -73°. 
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Fig. 8. — Comparison between 1 mm dust continuum data {upper row, as shown in Figure 
12]), model {middle row), and residuals {lower row) of L1157-mm in the image space. Images 
in each column share the same u-v imaging weighting and contour levels. The contour levels, 
noise rms (a), and beams are: column 1: [-3,3,4,5,7,10,14,20,30,42] xcr, a =4.0 mJy beam~^, 
1.77"xl.61" at a position angle of -38°; column 2: [-3,3,4,5,7,10,14,18,22,26] xa, a =5.5 mJy 
beam-\ 1.14"x0.98" at a position angle of -48°; column 3: [-3,3,4,5,7,10,13,16] xct, a = 7.0 
mJy beam~\ 0.71"x0.62" at a position angle of -51°; column 4- [-3,3,4,5,6,7,8] xa, a = 12.0 
mJy beam~^, 0.37"x0.30" at a position angle of -83°. The 3 a residuals at the smallest-scale 
image {lower-right) is likely due to the differences of the emission peak position measured 
using 1 mm and 3 mm data, while the modeling is done around the peak position measured 
using 3 mm data. If we shift the model with this offset at 1 mm, no residuals is left at 3 cr 
level in the small-scale image. 
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a are left in the small-scale image. 



The posterior- weighted results suggest a density power- law index p ~ 2. In this case, 
the envelope density structure is similar to a singular isothermal sphere or the beginning 
stage of the Shu model with a very small infall region. In the Shu model, a free-fall-like p ~ 
1.5 profile is established quickly during the collapse process. If the Shu model is applied 
strictly, an extremely young age of ~10^ yrs is implied. This age is much younger than other 
age estimates . For example, a kin ematic age of ~15,000 yrs is suggested by the outflow 
observati ons (iBachiller et al. 120011 ). The results are consistent with the single- wavelength 
study of iLooney et al. in which a larger sample of Class YSOs are modeled and 

unphysical young ages are derived using the simple self-similar model. A steep density profile 
can be related to a finite mass rese rvoir, as the constraint f rom an outer boundary can steepen 
the density in the outer envelope (IVorobyov fc Basull2005l ). Another possibility is the change 
of dust grain properties across the envelope. If larger grains are present towards the inner 
envelope, their greater opacity can result in a steeper density profile being estimated. We 
do not investigate the radial variation of /3 in our model, as no apparent /3 variation is seen 
in Figure O 



4.2. Spherical Power-law Envelope with an Inner Unresolved Component 

Disk formation is a natural consequence of angular momentum conservation when a 
rotating envelope collapses. It is expected to happen early in the star formation process, 
approximately in the Class stage. While characterizing disks in Class YSOs, in particular 
their size and mass, is critical to reveal the mass accretion process, observing them is diffi- 
cult due to the dusty envelopes around them. Distinguishing the disk component from the 
envelope emission requires a good understanding of the envelope, measurin g the unresolved 



emission as the circumstellar disk component (e.g., iKeene fc Massonlll990l : I Chandler et al. 



19951 ). 



Although the pure power-law envelope can fit the observational data with statistical 
significance ( §4.ip . we add another parameter, a point source fiux density, to the model to 
represent any unresolved component in our interferometric observations of L1157-mm. For 
simplicity, we assume that the dust properties for the unresolved component are the same 
as that for the rest of the envelope. Physically, this point source fiux density is interpreted 
as an upper limit of the embedded disk component with a size smaller than the highest 
observational resolution of ~0.3" or 75 AU. 



Using the same technique as discussed in §4.11 we characterize the model parameters. 
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Figure [9] shows the result marginahzed probabihty distributions, and Table |3] lists the ex- 
pectation values and uncertainties. Visualization of the model-data comparison using the 
marginalized parameters is shown in Figure [TUl Figure [TT], and Figure [T2] in the visibility 
domain and the image domain. 

While a better fit with a smaller is achieved by the the power-law envelope model plus 
an unresolved component, the results are consistent with the results of the pure power-law 
envelope model presented in §4.11 The posterior-weighted mean of the density power-law 
index p is slightly smaller than that of the pure power-law envelope model, due to the 
contribution of the point source flux. The dust opacity spectral index (3 and density po 
are consequently affected. With the added complexity of the model, larger uncertainties 
for the model parameters are obtained. In particular, the flux density of the unresolved 
component is not well constrained as it is not a necessary parameter. Nonetheless, the 
density index p is still close to 2, so the inconsistency with the Shu model still exists ( §4.ip . 
The posterior-weighted mean of the unresolved flux is ~2% compared to the total flux 
measured by single-dish observations (iGueth et al.ll2003l ). and ~11% of the flux measured 
by CARMA. 

The flux density of the unresolved component can be converted to the upper limit of the 
embedded disk mass within th e framework o f the m odel. If we follow the empirical method 
of disk mass approxira ation in iLooney et al.l ( 120031 ) based on the disk modeling of HL Tau 



m iMundy et al.l ( 119961 ) . a disk of 0.05 Mq at a distance of 140 pc is used as the standard 
candle for 100 mJy emission at 2.7 mm. As a result, our marginalized model gives a disk 
mass of 4.1 Mjup- Alternatively, we can use a single-temperature optically thin source model 
to estimate the mass, that is, 

M = FJ^KMr) (9) 



Table 3. Power-law Envelope Model with an Unresolved Component 



Parameter 


Mean 


Radius of the 68% confidence interval 






(statistical noise only) 


(with flux uncertainty) 


(1) 


(2) 


(3) 


(4) 


dust opacity spectral index (3 


0.84813 


0.00038 


0.11 


density at 100 AU (in 10"^^ g cm"^) . 


7.8627'^ 


0.0053 


0.96 


density power-law index p 


1.95004 


0.00031 


0.02 


unresolved 1 mm flux density (in mJy) 


19.1835 


0.0032 


1.5 



^corresponding to a total envelope mass of 2.0662 Mq 
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Fig. 9. — Same as Figure[5]for the spherical power-law model with an umesolved component. 
The absolute flux uncertainty is included. Marginalized posterior probability distributions 
for all four parameters (dust opacity spectral index (3, density po at 100 AU in the unit of 
10~^® g cm~^, envelope density power- law index p, and point source flux density at 1 mm in 
the unit of mJy) are shown. In the histograms, the dashed vertical lines enclose 68% or 1 a 
confldence interval, with the expectation values and a listed in Table [31 The dark and light 
areas in the 2-D contour plots are the 68% and 95% confldence regions. 
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Fig. 10. — Same as Figure [6] but for the power-law envelope model plus an unresolved 
component. Error bars are statistical errors in the annuli-bins only. Observational flux 
density, averaged vectorially and binned in u-v annuli around the source center, are shown 
by circles for the 1 mm data and asterisks for the 3 mm data. The model fit with the 
marginalized parameters is shown by solid lines, which includes two components: a power- 
law envelope (broken lines) and an unresolved disk (dotted lines). 
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Fig. 11. — Same as Figure [7] but for the power-law envelope model with an unresolved 
component. Comparison between 3 mm dust continuum data {upper row), model {middle 
row), and residuals {lower row) of L1157-mm in the image space. 
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Fig. 12. — Same as Figure M but for the power-law envelope model with an unresolved 
component. Comparison between 1 mm dust continuum data {upper row), model {middle 
row), and residuals {lower row) of L1157-mm in the image space. 
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( iHildebrandl Il983l ) . Following iLooney et al.l ( 120001 ) with the assumptions of T = 60 K and 
fi;=0.1(z//1200 GHz) cm^ (dust+gas), the estimated disk mass is 3.6 Mj„p. Although 
these two methods of disk mass estimation give consistent results, the mass estimate is 
subject to the uncertainty of dust emissivity and temperature (see §5.31) . 



4.3. Rotating Collapse Model 



Gr avitational collapse o f an en velop e with uniform rotation has bee n studied in lUlrich 



(119761 ) ■ iGassen &: MoosmanI (Il98ll ). and iTerebey. Shu, fc GassenI (Il984j . hereafter the TSG 



model). The initial condition of the TSG model is a singular isothermal sphere, as in the 
Shu model. The non-zero angular momentum causes material to fall onto the midplane, 
following the streamline equation 

L = , (10) 

Tc 1 — COS 6/ COS 6*0 

where Vc is the centrifugal radius, 6 is the angle from the rotation axis, and 6o is the angle of 
the streamline at large r. A disk structure is expected inside the envelope with the density 
distribution 

M cos 6* ^^^2/ cos 6^ 2cos^^^0x-i 

47r(G'Mr3)V2'^ ^ '^^J ^'^^^^ r/r, ^ ' ^ ' 

We adopt the TSG model for the envelope fitting. The model parameters include: (a) 
the dust opacity spectral index (3 as in Eq. (|3]), (b) the dust density po at 100 AU, (c) the 
centrifugal radius Vc of the TSG model in Eq. (fTOj) and Eq. (fTT]) . and (d) a point source flux 
density (at 1 mm) to represent any unresolved component. The unresolved component is 
assumed to have the same dust properties as the envelope for scaling the 1 mm flux density 
to 3 mm. Besides that the TSG model implies an embedded disk, this point source flux 
density is required to obtain good fits with the TSG envelope; we were not able to obtain a 
fit with 90% confidence level with a zero unresolved fiux. 

The model parameters are investigated as in §4.11 and §4.21 Figure US] shows the 
marginalized probability distributions of model parameters, and Table H] lists the expec- 
tation values and uncertainties. Furthermore, visualization of the model-data comparison 
using the marginalized parameters is shown in both the visibility domain and image domain 
in Figure [HI Figure [13 and Figure [TH] 

The TSG model only fits the data with a bright point source component. In this case, 
the envelope contributes little flux towards the total dust continuum, and the unresolved disk 
component dominates the emission at both 1 mm and 3 mm, especially at long baselines 
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Fig. 13. — Same as Figure [5] but for the TSC model with an unresolved component. The 
absolute flux uncertainty is included. Marginalized posterior probability distributions for 
all four parameters (dust opacity spectral index (3, density po at 100 AU in the unit of 
10^^® g cm~^, centrifugal radius Vc, and point source flux density at 1 mm in the unit of 
mJy) are shown. In the histograms, the dashed vertical lines enclose 68% or 1 a confidence 
interval, with the expectation values and a listed in Table HI The dark and light areas in 
the 2-D contour plots are the 68% and 95% confidence regions. 



- 30 - 



o 
o 
o 



>^ 

E o 



(0 

c 

Q 
X 

3 



100 



(arcsec) 
10 



1mm 



® 



.A A 



5^ 



N 



N 



\ 



\ * * 



— -\ 

3mm \ 



_I I I I I I I l_ 



_I I I I I I I l_ 



_l I I I I I L 



10 100 
uv distance (kA) 

Fig. 14. — Same as Figure[TO]but for the TSC envelope model plus an unresolved component. 
The total model fit (solid lines) includes two components: an envelope (broken lines) and an 
unresolved disk (dotted lines). 
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Table 4. TSC Model with an Unresolved Component 



Parameter 


Mean 


Radius of the 68% 


confidence interval 






(statistical noise only) 


(with flux uncertainty) 


(1) 


(2) 


(3) 


(4) 


dust opacity spectral index j3 


0.96008 


0.00022 


0.10 


density at 100 AU (in IQ-^^ g cm"^) _ 


2.7203 


0.0022 


0.34 


centrifugal radius Vc (in AU) 


44.9236 


0.0016 


14 


unresolved 1 mm flux density (in mJy) 


135.0764 


0.0067 


8.9 



'^corresponding to a total envelope mass of 5.0266 Mq 
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Fig. 15. — Same as Figure [7] but for the TSC model with an unresolved component. Com- 
parison between 3 mm dust continuum data {upper row), model {middle row), and residuals 
{lower row) of L1157-mm in the image space. 
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Fig. 16. — Same as Figure [H] but for the TSC model with an unresolved component. Com- 
parison between 1 mm dust continuum data {upper row), model {middle row), and residuals 
{lower row) of LI 157- mm in the image space. 
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(see Figure [M]) . While a rather shallow envelope model such as TSC is assumed, a strong 
small-scale c omponent is required to fit the data. S imilar results have also been seen in other 
studies (e.g., Terebey et al.lll993 : Enoch et al. As discussed in §4.2^ the unresolved fiux 



density can be converted to the upper limit of disk mass. The posterior-weighted parameter 
implies an embedded disk of ~25 Mjup using either method in §4.21 

Although the best-fit model passes a chi-square hypothesis test (or the null hypothesis 
is rejected with 90% confidence), the TSC model does not fit the data as well as the power- 
law envelope model does. As seen in Figure [IS] and Figure [111 the model does not subtract 
the data as cleanly as the power-law model does and leaves more residuals. In particular, 
residuals at 5 o" level are seen in the 3 mm image (column 2 in Figure [T5|) . A worse fit is 
also shown by the larger x^, which will be examined in more detail in §5.11 



Discussion 



5.1. Model Comparison 

A model is just a simplification of the unknown reality, but we want to know which 
model provides a better approximation to all available data. In this section we apply model 
selection techniques and rank the models. 

By applying Bayesian inference at the model level, one can use the Bayesian evidence 
for model selection. As discussed in Appendix B (Eq. f IBip and Eq. f lB2p ). the Bayesian 
evidence represents the probability of data given the model. It is marginalized over the 
full parameter space so the values of model parameters are not important, as opposed to 
parameter estimation for a particular model (Q. For comparing two competing models, the 
ratio of evidence, also known as the Bayes factor, represents p osterior odds and can infer 
whether one model is preferred over the other (see lLiddlell2009l . for a review). 



However, an exact method to compute the Bayesian evidence needs to fully evaluate 
likelihood in the entire parameter space and is very computationally expensive. The poste- 
rior distribution sampled by MCMC (§!]) peaks around the maximum posterior probability, 
and is not sufficient to calculate the Bayesian evidence. Approximation such as the use 
of informati on-theoretic methods is a good alternative approach for model selection (e.g., 
Liddlel l2007l ). For example, the Akaike information criterion (AIC, lAkaikd Il974j ). derived 
using the KuUback-Leibler information (or K-L distance), is defined as 

AIC = -21n£„„, + 2A;, (12) 

where C^ax is the maximum likelihood and k is the number of model parameters. AIC 
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provides a simple measure of how good the model approximates the information contained 
by the data, and a smaller value implies less information is l ost and hence a better mode l. 
Detailed derivation and statistical implications can be found in lBurnham fc Anderson! (120021 ). 
A second-order AIC, or AIC corrected, is suggested for small-sample bias adjustment as in 



AIC. = AIC 



2k{k + l) 
N-k-r 



(13) 



where N is the number of data points. But in our case, ^ k so AIC and AICc converge. 
On the other hand, the Bayesian information criterion (BIC. ISchwarzlll978h . defined as 



BIC = -2lnC^ax + k\nN, 



(14) 



is an approximation based on the Bayesian evidence ratio (also see iLiddld l2004l ) . 



A good model seeks for balance between goodness of fit and model simplicity. To obtain 
a better fit to the data or a smaller x^, one may increase the model complexity with more 
parameters, but unnecessary use of parameters and over-fitting should be discouraged. The 
tradeoff is also seen in Eq. (1121) and Eq. (I14p . As the best model minimizes AIC and BIC, 
smaller decreases the first term but extra parameters increase the second term. Model 
complexity in terms of the number of parameters is penalized in either AIC or BIC. 

In this study, we evaluate AICc and BIC for all models. Because either AICc or 
BIC is on a relative scale, only the differences instead of actual values are meaningful 
(IBurnham fc Anderson! |2002| ). Results are listed in Table !5l, where the model with the 
smallest value is preferred. As AICc and BIC suggest different ranking between the pure 
power- law envelope model and the power-law envelope plus an unresolved component, we 
do not make an inference between them. However, the results show a large positive AAICc 
and ABIC for the TSC model, implying that the power-law envelope model (either with or 
without the unresolved component) is decisively preferred against the TSC model plus an 
unresolved component. 



Table 5. Model Comparison 





Model 


# of parameters 




AAICc 


ABIC 




(1) 


(2) 


(3) 


(4) 


(5) 


Power-law envelope 




3 


9613473.1 


1.9 





Power-law envelope 


-|- an unresolved component 


4 


9613469.2 





12.2 


TSC envelope -|- an 


unresolved component 


4 


9615698.7 


2229.5 


2241.7 
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The model selection results are consistent with the a posteriori check presented in §H 
Compared to either the pure power-law envelope model or the power-law envelope plus an 
unresolved component, considerable residuals are seen in the image domain for the TSC 
model plus an unresolved component. 



5.2. Grain Growth 



As mentioned in §3.3[ dust properties are characterized by the opacity spectral index /3 
in our modeling. Depending on the environment, /3 typically varies between and 2. While 
/3 ~ 2 implies sm all grains as in th e interstellar medium, a smaller /3 is usually found in 
many YSOs (e.g., iNatta et al.l 120071 ). Decrease of (3 can be caused by many factors, such 
as change of composit ion or grain geometry, but is u sually associated with change of grain 
size distribution (e.g., iKriigel fc SiebenmorgenI 119941 ). The (3 value can be an evolutionary 
indicator of the dust grains in YSOs, as small grains in YSO s grow through coagu l ation , 
and eventually form planets if conditions allow. Nevertheless, iMiyake fc Nakagawal (119931 ) 
studied the size effect and showed that the observed decrease of /3 in disk regions can be 
explained by the growth of grain size without change of chemical composition. If the grains in 
a dense region are composed of the same materials as in the interstellar grains, the maximum 
grain size is e xpected to be larger than 3 mm to explain the observational results of /3 < 1 
( iDraind l2006l ). Grain growth has bee n the most widely ac c epted explan ation for the small 
/3 found in protoplanetary disks (e.g., iBeckwith et al.ll2000t iDraind 120061 ). 



In the three envelope models presented in §U the posterior- weighted mean of /3 range 
from 0.84 to 0.96. The absolute flux calibration dominates the uncertainty of /3 estimation, 
resulting in a systematic uncertainty of ~0.1. Still, (3 being significantly smaller than the 
interstellar value is indicated. Our estimate fo r L1157-mm i s in a greem ent with the 



sample s of Class YSOs in |j0rgensen et al.l (120071 ). iKwon et al.l ( 120091 ) . and IShirley et al. 



(l2011bl ). The result implies that dust grains in LI 157- mm, and arguably most Class YSOs, 
have gone through some grain growth to at least millimeter size from the initial interstellar 
grains. However, when ex actly dust grains st art to grow during the protostellar evolution 
is uncertain. For example, iRicci et al.l ( l2010al ) compared /3 of YSOs with their evolutionary 
ages and did not find apparent trend or difference in /3 for YSOs in different evolutionary 
stages. 

As a uniform dust grain property is adopted in our simple dust model, our estimate of (3 
represents the grain property in the whole system, including disk and envelope. Depending 
on the assumed envelope model, the embedded disk can contribute a significant fraction of 
flux. Compared to grains in the envelope, grains in the disk are expected to be larger in 
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size as an initial step of planet formation. Therefore, a smaller /3 is expected in the disk 
than in the envelope. Besides, grain properties are likely to var y across the envelope, as 
has been observat ionally suggested for some Class sources (e.g., IChandler fc Richei]l2000l: 



Kwon et al.l |2009[ ) . Since our L1157-mm data are consistent with a single /3 value accross 



the envelope, the radial dependence of /3 or a disk with a different /3 is not modeled in this 
study; to address the dust property change in YSOs requires a more complex model. 



5.3. The Earliest Circumstellar Disks 

Circumstellar disks form as a physical consequence of angular momentum conservation 
when protostars accrete materials from their surrounding envelopes. These planet-forming 
pr e-main-sequence disks h ave been observed and studied extensively (e.g., see the reviews 
of IWilliams fc Ciezall201ll ). In particular, disk evolution from early Class to Class I stage 
is interesting as it is the phase during which most mass accretion occurs. The mass and 
size of the disk grow rapidly as the system evolves and depend on the rotation rate and the 
magnetic field strength of the backg round cloud, as well as the detai l ed mechanisms of angu- 



lar n i qmentum redistri bution (e.g., iTerebey et al.lll984j : iBasulllQQSt iMachida &: Matsumoto 



201ll : lDapp et al.ll2012[ ). 



While the mass and size of these youngest disks are essential to reveal the early pro- 
cess of disk formation, observing them is, however, not straightforward. In addition to 
the limitations of observational resolution and sensitivity, these disks are deeply embedded 
in their natal envelopes; probing them usually relies on indirect methods. For example, 
detection of water and methanol lines in Class protostars can constrain the embedded cir- 



envelope interface, while different scenarios may not be ruled out (e.g.. 


Goldsmith et al. 


1999; 


Velusamv et al. 


2002; 


Watson et al. 


2007; 


J0reensen & van Dishoeck 


2010). Near- infrared 



scattered light images, showing a dark lane along th e edge-on disk, have also been used to 
infer the embedded disk structure (ITobin et al.ll2010l ). 



Another way to probe these youngest disks in embedded YSOs is through observa- 
tions of dust continuum with a two-component model. The millimeter data consists of 
data continuum from both the disk and the envelope; with accurate modeling of the enve- 
lope, the circumstellar disk comopnent can be separated. In other words, the disk com- 
ponei it is measured as the residual emission with the envelope contribution subtracted 



Keene &: Masson 



Chiang et al. 



2008 



1990 



Lqonev et al.l l2003l : iHarvey et al.l l2003l : jjqrgensen et al.l 12005 



Enoch et al.ll2009l ). This method avoids the complexity of chemical ef- 



fects, but requires the use of a theoretical envelope model. 
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Using a two-component model, attributing the entire unreso l ved fl ux t o be from the 
embed ded disk, and following the mass estimation in iLooney et al.l ( 120001 ) and iLooney et al. 



( 120031 ). we derive the disk mass to be ~4 Mj^p assuming a power-law envelope, or ~25 M 



Jup 

assuming a TSC envelope. Note that this mass estimate is an upper limit valid only within 
the framework of the models, and is highly dependent of the assumed envelope structure, dust 
opacity, disk temperature, and optical depth. For the envelope structure, we have shown that 
the power-law envelope model is preferred against the TSC model. But a large uncertainty 
still exists. For example, if instead we adopt a different dust opaci ty ft:=0.015(i//3C )0 GHz ) 



cm^ g ^ (dust+gas) and a single temperature of 30 K, as used in iGreaves &: Ricd ( 120111 ) 



the disk mass is ~13 Mj^p assumin g the power-law envelo pe model. Our deduced disk mass 
is low compared to the estimate in iGreaves &: Ricd (1201 11). where a disk mas s of 80 M j„p is 
obtained for L1157. Although Eq. is also used in IGreaves &: Rice ( 2011 ). they assume 
the compact emission at 2" is solely from the disk and no envelope subtraction is done, which 
possibly causes the differences in the estimated disk mass. 

On the other hand, an empirical method that measures the small-scale flux at baseline 
~50 kA has also been used to estirnate th e unresolved disk component in embedded YSOs 
( 1j0rgensen et al.ll2009t lEnoch et al.l 1201 ll ). This is based on the presumption that the en- 
velope contributes little flux >50 kA. The 1 mm flux density of L1157-mm i s around 200 
mJy at 50 kA, implying a disk mass of ~174 Mjup using Ki.3mm=0.009 cr n^ g~^ (lEnoch et a 



(120111 ) 



201ll ). The estimated dis k mass is comparable to other Class sources in lEnoch et al ^ , 

and lighter than those in I J0rgensen et al. J2009h . but much heavier than our disk mass es- 
timate using a two -component model. The empirical method of |j0rgensen et al.l (120091 ) and 
Enoch et al.l (12011! ) seems to over-estimate the disk component in L1157-mm; one reason is 
that the flux density drops signiflcantly longward of 50 kA (Figure [1]) so the flux from the 
unresolved disk is apparently lower than 200 mJy. In ad dition, a relatively shal low envelope 
proflle is assumed in either study (power-law of 1.5 in |j0rgensen et al.l l2009l and TSC in 
Enoch et al.l 1201 ll ). and the assumed envelope structure considerably affect the flux ratio 
from envelope and disk. 



Since our observations do not resolve the circumstellar disk, we can put an upper limit 
on the disk size for L1157-mm. Assuming a distance of 250 pc, the u pper limit o f the d isk 
radius is ~40 AU. This is consistent with the theoretical scenario of iDapp et al.l (120121 ) in 
which Class disks are small. This size constraint for L1157-mm is also consistent with 
those in o ther Class YSOs. For example, no disks are detected for a sample of 5 Class 
objects in iMaury et al.l (120101 ) with a resolution down to 0.3", corresponding to a size scale 
of 42-75 AU, and the embedded dis k at the edge-on Class YSO VLA 1623A is constrained 



to be smaller than 50 AU in radius (I Ward- Thompson et al.ll201ll ). On the contrary, the disk 



embedded in the edge-on Class YSO L1527 has been resolved by 7 mm VLA observations 
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m 



Loinard et al.l ( 120021 ): the disk structure is also seen with SMA and CARMA observations 



(Tobin et al. 2012). The size of Class disks is comparable to or sma l ler than the size 



of old er circumstellar disks (e.g., lEisner et al.l 120051 : [Andrews et al.l l2009l : IVicente &: Alves 



20051 ). but no clear trend can be inferred at this point. 



6. Summary 

1. Multi-configuration CARMA observations of the edge-on Class YSO L1157-mm are 
presented. In our dust continuum data at both 1 mm and 3 mm, a nearly spherical 
circumstellar envelope is seen at the size scale of ~10^ to ~10^ AU. No circumstellar 
disk on the small scale is resolved. 

2. Radiative transfer modeling is performed to compare the interferometric data with the 
theoretical envelope models. A power-law envelope and a TSC envelope are considered. 
We add an unresolved component to represent the embedded disk. Bayesian inference 
is employed for parameter estimation. The absolute amplitude uncertainty, resulting 
from the flux calibration of the data reduction process, plays a critical role in parameter 
errors. 

3. A density index p ~ 2 i s sug gested for the power- law envelope, consistent with the 



results in iLooney et al.l (120031 ) for a larger sample of Class YSOs. An unphysical 
young age is suggested if the Shu model is applied strictly. The data can be fitted 
by a pure power-law envelope without a compact emission from the embedded disk 
component. 

The dust grain properties of the envelope are studied through the dust opacity spectral 
index /3. The result /3 ~ 0.9 is significantly smaller than the /3 value in the interstellar 
medium, implying that grain growth has already started in LI 157- mm. 

The unresolved disk component is constrained to be <40 AU in radius and ^4-25 Mjup 
in mass. However, the mass estimate of the embedded disk component heavily relies on 
the assumed envelope model as well as the assumed disk characteristics. For example, 
a shallow envelope, such as the TSC model with a density power-law index p ~ 1.5 in 
the outer region, requires a strong point source flux from the unresolved disk, while a 
steep envelope with p ~ 2 can fit the observational data without an embedded disk. 

Different envelope models are compared using an information-theoretic approach. The 
results prefer the power-law envelope model against the TSC model, which is also 
shown in the a posteriori check in the image domain. 
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7. This is the first study that utihzes the Bayesian techniques and model selection to 
consider multiple envelope models and make statistical inference for embedded YSOs. 
Future observations, especially high-resolution ALMA observations, will resolve the 
transition zone between the envelope and the disk, and further constrain the structures 
of Class YSOs. 
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A. Error Estimate of the Approximate Dust Opacity Spectral Index 



In the optically thin limit and the Rayleigh- Jeans regime, the dust opacity spectral 
index Ptun can be approximated using the fiux density at two wavelengths. In this appendix 
we discuss the error propagation from the observational uncertainty to the deduced PtMn 
value. Let Fi and F2 be the fiux density at frequencies vi and z/2, Pthm can be expressed as 
in Eq. 

A,.„ = Y'-'f' - 2, (Al) 

In Ui — In 1^2 

Assuming Fi and F2 are indenpendent variables with standard deviations cfi and (T2, the 
standard error propagation gives 



a 



l^thin 



dPthin 


to 


dPthin 


dFi 




dF2 



(To 



(A2) 
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Taking the partial derivative of Eq. flAl|) . we obtain 

dPthin 1 



dFi (In vi — In 1^2) Fi 

and 



(A3) 



dF2 (Inz/i -lnz/2)F2' ^ ' 

Replacing Eq. (IA2p using Eq. (IA3p and Eq. f lA4p . the uncertainty of the derived Pthin is then 

2 / 1 \ / ^1 , 0^2 



" [ln,.,-\nu2 ) [f^^fIJ- ^^^^ 

Using Eq. (lASP and assuming the absolute flux error can be represented as a Gaussian 
noise with a standard deviation of 10%, that is, ai = O.lFi and o"2 = O.IF2, the uncertainty 
of Pthin is ~0.15 for our data at 229 and 91 GHz. Note that here the error is assumed be 
normally distributed, different from the flat prior assumed in §4. 



B. Fitting Technique and Statistical Inference 

In this appendix we give a brief introduction to Bayesian inference, as opposed to 
frequentist statistics. Also, we describe our technical procedure to characterize model pa- 
rameters and their uncertainty. 

The main concept of Bayesian inference is to incorporate prior knowledge on the hy- 
pothesis. Also, information is represented in terms of a probability density function (PDF) 
in parameter space. Mathematically, given the observed data, the posterior probability of 
model parameters can be specified by the Bayes' theorem 

where x stands for model parameters, D stands for data, M denotes a particular model with 
its model assumptions and other background information, P{x\M) is the prior probability 
of model parameters x, P{D\x^ M) is the conditional probability or the likelihood of data 
given the model with parameter x, and P{D\M) is the evidence or global likelihood. The 
evidence P{D\M) is the net probability of the data given the model, as it sums the product 
of likelihood and prior over parameter space: 

P{D\M)= j P{D\x,M)P{x\M)dx. (B2) 
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The evidence is independent of the parameter values and can be seen as a normahzing factor 
in Eq. f lBip : therefore it is not important for parameter estimation of a single model. For the 
same reason, the model label M is sometimes omitted when only one model is considered. 
However, evidence is useful for comparing multiple models ( §5.11) . 

Whether to view a statistics problem with Bayesian or frequentist approach is under 
de bate. The disputes are beyond the scope of this study and more discussions can be found 



m 



Loredd (jl990l . Il992l ). Nevertheless, in the case of a uniform prior, the method of using the 
posterior probability is equivalent to maximum likelihood estimate as far as identifying the 
best-fit parameter values is concerned. As for full parameter estimation including estimating 
the parameter errors, different approaches are adopted for Bayesians and frequentists, and 
will be discussed later. 

The likelihood of data given the model is characterized by (as defined in Eq. ([7])) and 
P{D\x) ~ exp(— x^(x, D)/2) with model parameter x and observational data D. Therefore 
the most probable parameters with the maximum likelihood can be obtained by locating 
the global minimum of x^. We use the Nelder-Mead simplex algorithm as implemented in 
MATLAB with bound constraints to search for the minimum. Besides fast convergence, 
this method does not evaluate function derivative, which suits our application because fewer 
modeling evaluations are required. Several starting points are used to look for several con- 
vergent minimums, and they are checked to be consistent with each other. This is to make 
sure that what is found is the global minimum, not a local minimum. 

Once the best-fit parameter values are identified, efforts are made to characterize the 
uncertainty. The essence of parameter estimation is to characterize the reliability of an 
estimate on model parameters under the assumption that the best-fit model is correct. In 
other words, the best-fit model needs to show statistical significance based on a hypothesis 
test before any of the following parameter estimation can make sense. For example, in the 
standard Pearson's chi-square hypothesis test, value of the model needs to be smaller 
than a critical value depending on the degrees of freedom to reject the null hypothesis. In 
the following we discuss two common methods of parameter estimation: (1) a frequentist 
approach to characterize Ax^ and infer statistical significance, and (2) Markov chain Monte 
Carlo (MCMC) in the context of Bayesian inference. 

The frequentist Ax^ statistics has been suggested in Lampton et al. (1976) and Avni 
(1976), and summarized in Press et al. (2002). With the definition of Ax^{x) = x'^i^) ~ 
X^{^ best- fit), Ax^(x) is chi-square distributed with p degrees of freedom, where p is the 
number of fitted parameters or parameters of interest. Then the level of confidence can be 
estimated according to the chi-square distribution. Although it relies on the validity of the 
best-fit model, the exact value of x'^i^best-fit) is not important for parameter estimation. 
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The Ax^{x) statistics is independent of the Pearson's chi-square test and x^i^) statistics, 
and focuses on the variation of x^i^) i^i parameter space x. A common way to illustrate the 
results is through iso-chi-square contours or hyper-surface in mult i- dimensional parameter 
space as the confidence region. In the case that only partial parameters are of interest, 
the remaining nuisance parameters should be varied to minimize Ax^ix) instead of direct 
projections, (c.f. In Bayesian inference, the nuisance parameters are marginalized over.) 
For example, when only one parameter is of interest, Ax^(x) is distributed as a chi-square 
distribution with one degree of freedom. The 68% confidence interval corresponds to the 
region bounded by Ax'^{x) = 1. 

Despite t he controversy over the flaws of applying this method with nonlinear models 
(ILoredol Il992l ). estimating Ax^(x) over a large parameter space can be computationally 



difficult. A grid on parameters or equivalent technique is required. The large number of 
evaluati ons usually makes this method impractical, especially when the number of parameters 



is large (lFordll2005l ). 



On the other hand, MCMC offers a very efficient way to estimate the posterior prob- 
ability in Bayesian inference, compared to any other methods that require grid searching. 
Rather than minimizing on each grid point and probing the variation of x'^i^)^ the posterior 
probability respect to parameters of interest is estimated through marginalization over all 
other parameters. For example, given a PDF P{xi, X2I-D) where xi is the parameter of inter- 
est and X2 is a nuisance parameter, the X2 space is integrated over according to probability 
to obtain the marginalized PDF, as in 



P(x 



l\D) = J P{xi,X2\D)dX2 = j P{Xi\x2,D)P{x2\D)dX2. (B3) 

At first glance, a straightforward marginalization can be very computationally extensive, 
similar to the necessity of a grid evaluation in frequentist methods. However, marginalized 
results can be obtained efficiently with MCMC. One of the reasons is that it searches the 
parameter space according to probability and the parameter space with low probability is 
less explored and sometimes not probed at all. 

The Metropolis-Hastings algorithm of the MCMC method is utilized to construct the 
Markov chain. Markov chain is a sequence of parameter values representing the system 
and characterized by a transition probability that controls the random process from one 
state to another. The transition probability and the next state are only dependent of the 
current state, but not any previous states. Regardless of the starting state, the chain even- 
tually converges to a stationary or equilibrium distribution according the PDF. We use the 
Metropolis-Hastings algorithm to draw the sample and construct the chain. This algorithm 
uses a proposal distribution g(x'|x), or the candidate transition probability distribution func- 
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tion, to generate a trial state x' based on the current state x. Then the proposed state is 
randomly accepted with the acceptance probability 

' P{x'\D)q{x\x') 



a[x \x] 



mm 



P{x\D)q{x'\x) ' 

or otherwise rejected. The arrangement results in a transition probability 

T{x'\x) = q{x'\x)a{x'\x) 



(B4) 



(B5) 



which is reversible {7i{x)T{x'\x) = 7r(a;')T(a;|a;'), where tt{x) is the equilibrium probability at 
state x) and irreducible (possible to go from any state to any state). As introduced earlier, 
P{x\D) is the posterior PDF given the observational data, and approximately proportional 
to exp(— x^(a;, D)/2) with flat prior. Practically, x^i^^ D) is evaluated at each proposed state 
change. 

The Metropolis-Hastings algorithm assures that the chain converges to P{x\D) as the 
sample number is large. The convergence rate is related to the choice of q{x'\x). A typical 
choice is a Gaussian function centered around x, that is. 



q{x'\x) 



exp 



ix' — xY 



2w^ 



q{x\x'). 



(B6) 



The width of the Gaussian, specified by w, determines the trial step size. If the step size is 
too large, most trial states are rejected so the calculation becomes very inefficient; if the step 
size is too small, the chain behaves like a random walk and requires a long time to converge. 
An opt imal acceptance rat e is suggested to be around 0.23 for mult i- dimensional parameter 
space (iGelman et al.l 120041 ). The choice of a symmetric proposal distribution also reduces 
the acceptance probability into a simper form 



a[x \xi 



mm 



P{x'\D) 
_P{x\D) 



mm 



exp 



x\x,D)-x\x',D) 
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The posterior probability is obtained with a converged Markov chain, as its density of 
points in parameter space follows the posterior probability of the parameters. Marginaliza- 
tion is done through projecting the Markov chain to the space of parameters of interest. 
We estimate the 68% and 95% confidence limits and the corresponding standard deviation 
based on the simulated MCMC. Specifically, the 68% or 1 a confidence hmit encloses 68% 
of accepted points along the Markov chain, and represents the region containing 68% of 
the total probability distribution. We also report the expectation value of each parameter, 
weighted by the posterior marginalized probability as in 

(xi) = / P{xi\D)xidxi = -^T^jXij (B8) 
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where j denotes the points in the Markov chain and N is the total number of points. The 
expectation values from the marginalized distributions do not need to be identical to the 
parameters with the maximum likelihood, because we are projecting the values from a high 
dimensional distribution which may not be a multivariate Gaussian; however, they should 
be consistent. 
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